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Abstract 



We solve the generalized Langevin equation driven by a stochastic force with power- 
law autocorrelation function. A stationary Markov process has been applied as a 
model of the noise. However, the resulting velocity variance does not stabilizes 
but diminishes with time. It is shown that algebraic distributions can induce such 
non-stationary affects. Results are compared to those obtained with a deterministic 
random force. Consequences for the diffusion process are also discussed. 
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Modeling a physical system in terms of Langevin formalism must take into account the 
nature and origin of the stochastic force. Usually that force is taken in the form of the white 
noise but in many cases that is an unrealistic idealization. Among the systems possessing 
a finite noise correlation time, those with power-law (algebraic) correlations are especially 
interesting because of lack of characteristic time scale and divergent moments. Such systems 
are not unusual. The algebraic random force autocorrelation function (FAF) appears in the 
fluid dynamics |IJ,[J and linearized hydrodynamics [EJ. For such phenomena as the noise- 
induced Stark broadening [[| and nuclear collisions || , correlation functions proportional to 
1/t have been found. The latest form of correlations is of special importance for molecular 
dynamics because it corresponds to the problem of scattering inside a periodic lattice || . 

For systems with finite noise correlation time, the ordinary Langevin equation must be 
generalized [pfl to ensure proper fluctuation-dissipation relations. In the absence of external 
potential, this equation has the form 

m MQ_ = _ m r l K ^ _ T ^ dT + F{€) (1) 



dt Jo 

where F(t) is a stochastic force and K(t) denotes the retarded friction kernel. The 
fluctuation-dissipation theorem imposes the relation ||: K(t) = (F(0)F(t)) /mT, with tem- 
perature T and mass m. We require (F(t)) = and the following FAF: 

O,(t) S <F(0)F(«)>={^ £ ^ I (2) 

where e is a small number. The coefficient (3 we set equal to one. The solution of Eq.([TJ) 
with the initial condition v(0) = takes the form [0] 



v(t)=g(t)+ [ R(t-r)g(r)dr (3) 

where g(t) = m~ l f F '(r) dr and the resolvent is given by R(t) = exp(— at) (c\ sin bt + 
C2Cos6t) +m / °°xexp(— tx)/[(mx+Eix(ex) — l) 2 + 7r 2 ] dx. The modified integral exponential 
function Eii(x) is defined by the series: Eii(a;) = 7 + lnx + Y^Li x n /n\ n, where 7 = 
0.5772157. . . is the Euler constant. The other constants are fixed in the following at the 
values: a = —3.52832, c\ = —4.76673, c 2 = —5.35498 and b = 2.49975, corresponding to 
e = 0.01 and m — 1. 

We are interested in evaluation of the second moment of velocity distribution (v 2 ) (t) . It 
can be obtained directly from (0): 

(v 2 )(t) = 2&dT(t - t)C f {t) + 2^dr^d Sl j:ds 2 R(t - r)C F (\ Sl - s 2 \) 

+ J * dr J* ds % d Sl Jo S ds 2 R{t - T )R{t-s)C F {\ Sl -s 2 \). (4) 

Some of the above integrals have to be performed numerically. The result for T = 1 is 
presented in Fig.l. As expected, the system reaches the equilibrium state. 

In the following, we will refer to the above result as "analytical". Alternatively, we 
can introduce a concrete stochastic process characterized by covariance in the form (Q) and 
calculate (v 2 ) from a Monte Carlo simulation. We apply the "kangaroo process" (KP). 



It is defined |10| as a stepwise random function: F(t) = Fi = const in the time interval 
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ti <t < ti+\. The length of interval of constant F, s, is a function of the value of the process 
itself. The KP is a stationary Markov process, determined by a stationary probability 
distribution Pkp{F). It can be easily defined for arbitrary covariance. We get [O the 
required form (fj) by choosing P KP (F) = l/(2a) = const for F G (—a, a), zero elsewhere, 



where a = yS/rj with rj = e/T. The time increment corresponding to a given F follows 
from the formula s = 3a\F\~ 3 . Thus s assumes values between 77 and infinity and its density 
distribution is of the form 

P(s) = (r ] 1 ' 3 /3) s- 4 / 3 6(s-ti), (5) 

where 6{x) is the step function. Moments of F can easily be obtained by averaging over 
the uniform distribution P KP (F). We get (F) = and Var (t) = (F 2 ) = l/rj. Since 
2Pkp{\F\)(1\F\ = P(s)ds, we can average also over P(s): 

roc 

Var (t) = 3r]- 1/3 / S - 2/3 P(s)ds. (6) 

Jr] 

Also the KP covariance |T0| , |TTH Ckp can be expressed in terms of the distribution P(s) 



C KP {t) = 3?T 1/3 / s-^ 3 exp(-t/s)P(s)ds. (7) 

Jr) 

Inserting F(t) into Eq.(|3D allows us to determine the velocity time series of the Brownian 
particle. The variance at a given time t is obtained simply by calculating v(t), squaring it 
and averaging over many trajectories. Fig.l presents the result: (v 2 ) does not stabilizes at 
the expected value (v 2 ) = T/m but instead it dwindles with time, obeying the approximate 
relation (v 2 )(t) ~ t~ 0S7 . 

The above outcome is surprising because the stationary process brings about an appar- 
ently non-stationary result. Moreover, according to (|J), (v 2 ) is completely determined by 
the covariance of F and every simulation satisfying @) should reproduce the analytical re- 
sult. In order to understand the origin of that inconsistency, let us reconsider in details how 
actually the stochastic force value enters the Langevin equation. Evaluation of the Brownian 
particle velocity requires the value of F at a given time t. For that purpose the distribution 
of s is crucial because this value follows from the length of current interval in the stepwise 
evolution of KP. However, the requirement that we choose only those intervals which con- 
tain the time t imposes some bias, e.g. longer intervals are more probable. Therefore a 
distribution we actually use in the Langevin equation (the "effective" interval distribution), 
may not be identical with P(s). Its cumulative distribution function, $(s,t), can be derived 
in the following way. First let us consider s < t and assume that t is found in n + 1 interval, 
i.e. S n = si + S2 + . . . + s n < t and 5* n+ i > t. The probability that the sum of n intervals 
yields a value between x and x + dx we denote by P n (x)dx, providing that each component 
has the distribution P(s). The distribution function <&(s, £) is just equal to the conditional 
probability that an interval is larger than t — x, for any x between t — s and t, and any n from 
1 to N, where iV denotes the integer part of t/r]: = Yln=x Its S n (x)dx J^_ x P(£)d£. 

Introducing S(x) = Y<n=i Sn{x) an d inserting P(x) from (||), we get the following equation: 

$(s,t) =r] 1/3 S(x)[(t-x)- 1/3 - s~ 1/3 ]dx for 7]<s<t. (8) 

Jt-s 
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For s > t the lower limit of integration extends to zero. Moreover, we have to take into 
account also events for which t is contained already in the first interval. The final formula 
reads: 

$(s,t) =V /3 ( t S(x)[(t-x)- 1/3 -s- 1/3 ]dx + r 1/3 -s- 1/3 \ for 8 >t. (9) 



The direct evaluation of S(x) is very difficult. We can avoid it utilizing the normalization 
condition $(oo, t) = 1. The function S(x) must then satisfy the integral equation 

S(x)(t - x)- 1/3 dx + r 1 ' 3 = ?T 1/3 , (10) 
called Abel's equation. It possesses a weakly singular kernel, depending only on the difference 
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of its arguments. Therefore we can apply the Laplace transforms technique to solve it [|12 
The solution is of the form 

S( x ) = c r r]- 1/3 x~ 2/3 - 5(x), (11) 

where a constant cr = l/[r(l/3)r(2/3)] ~ 0.2757... contains the Gamma function. In- 
serting S(x) to (|S]) and ([]) and evaluating integrals gives us the expression for <E>(s, £). To 
obtain the probability distribution V(s,t), we have to differentiate $(s,t) over s. The final 
result is simple: 

( , _f Cr [t 1 /3_ (t _ s) 1/3] ,-4/3 for ^< s < t 

{ , > \(c r * 1 /3 + ^1/3/3) s -4/3 f or s>t . 

The effective interval distribution appears to be time-dependent and consisting of two 
branches which do not join smoothly. We encounter a similar problem asking about the 
mean time we must wait for a bus, providing we know the average time interval between 
subsequent bus arrivals (t). The answer is not r/2, as one could expect, but just r. This 
"waiting-time paradox" [13[] can be elucidated by calculating the effective, time-dependent 
probability distribution, analogous to (12). In that case, however, the original distribution 



is an exponential which results in the fast equilibrization and the effective distribution 
asymptotically becomes time- independent. For V(s,t) it never happens. Moreover, since 
the probability V(s > t,t) = / t °° V(s,t)ds = 3cr ~ 0.83 does not diminish with time but 
remains constant, the entire distribution shifts with time towards long intervals. In fact, 
this outcome is not unexpected because all moments of V(s,t), as well as of P(s), are 
divergent. On the other hand, long intervals correspond to small values of the process itself, 
which points out a reason of declining of the variance. To derive expression for the effective 
variance V(t), we can use Eq.(|]) substituting V(s,t) for P(s). Evaluation of the integral 
gives us the final formula 

V(t) = c r rT 1/3 [3 In 3/2 + nV3/6 + ln(t/r})\ r 2/3 (t^rj). (13) 



Hence the variance really decreases with time [14]. The KP covariance must also be modified. 



Replacing P(s) in Eq.(0) by V(s,to), where to is an initial time, and evaluating integrals we 
get the effective covariance 
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C K p(t, t ) = cr^Hl'h- 1 [3 - exp(-t/2t )r(l/3) W- 1/3 ,- 1/2 (t/t )] , (14) 

where W a ^ (x) is the Whittaker function ||15|| . This result is quite different from the orig- 
inal covariance (Q) and explains us why the simulation does not agree with the analytical 
prediction (|J). CKp(t,t ) ~ t~ l for large t but it depends also on to- 

The above conclusions imply that problems involving algebraic correlations (e.g. critical 
phenomena, hydrodynamics), investigated in the framework of the Langevin description, 
should be handled with caution. Conversely, an experimental evidence of declining vari- 
ance in such systems does not necessarily mean that Langevin formalism obeying standard 
fluctuation-dissipation theorems does not apply. In any individual case one should examine 
the distribution P(s), the shape of which for large s decides whether the system behaves in a 
stationary way. Is the stationary behaviour possible at all for correlations (|2])? The analyt- 
ical result would be valid for a steep P(s). The fastest decaying distribution one can obtain 
with the KP for (|2]) declines asymptotically like s~ 2 Since also for this distribution all 
moments diverge, we expect similar affects as for (|5|). 

Another possibility is to apply a deterministic process, instead of the Markovian stochas- 
tic one. For that purpose, let us consider a two-dimensional lattice of periodically distributed 
disks of radius r, with a particle bouncing elastically from them. Then the particle motion 
is free between subsequent collisions and its velocity u = (u x ,u y ) = const. This system, a 
periodic Lorentz gas, is equivalent to the Sinai billiard with periodic boundary conditions. 
We assume 2r < I, where I is the distance between disks centers. The system is strongly 
chaotic but the autocorrelation function of either component of particle velocity falls off 
slowly, as 1/t for large t |J. Therefore we can simulate solutions of (p]) assuming the ve- 
locity of particle inside the independently evolved Sinai billiard as the stochastic force F(t) 
fT6|| . A quantity of interest is the distribution of free paths: it falls like s~ 3 |T7j , steeper then 



for any KP. Its mean is convergent and the second moment weakly divergent. For numerical 
simulations we assume / = 1, r = 0.8, |u| = 1 and F = 7.3VTu x . Then Cp = T/t for large 
t. We must stress, however, that the form of FAF at small t also influences solutions of (P. 
Thus the simulations utilizing the Sinai billiard should be regarded as an approximation. 
The result of the numerical calculation of the variance {v 2 )(t) presents Fig.l. There are some 
discrepancies at small t, comparing to the analytical prediction, that can be attributed to 
differences in FAF. Asymptotically however, (v 2 ) stabilizes at the equilibrium value and 
both results coincide, in contrast to the KP case. 

Finally, we wish to calculate the velocity autocorrelation function (VAF) C v (t) = 
(v(to)v(to + t)), which is responsible for transport properties of the system. It allows 
us, namely, to determine the diffusion coefficient T> = J Q °° C v (t)dt. Typically, T> is finite 
which corresponds to the normal diffusion. The analytical result for VAF in our case is 
0: C v (t) = T/m [1 + Jq R{t)(It\. We present this function in Fig. 2. It has the tail of 
the power-law shape with numerically estimated exponent equal to —1.18. On the other 
hand, we have calculated VAF from simulation utilizing the KP ||18|| . The result for two 
values of to is presented in Fig. 2. We have normalized both functions to unity at t = 0. 
Their shape is very different from the analytical result. The VAF initially falls but then it 
stabilizes. It depends strongly on to, becoming more flat for larger to- The stationary case 
applying the Sinai billiard also produces result different from the analytical one - C v (t) is al- 
ways non-negative and does not approach zero for increasing time. As regards the transport 
properties of the system, the analytical C v (t) implies the normal diffusion. Determination of 
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the precise shape of VAF at large t for the KP and the simulation utilizing the deterministic 
random force, requires further studies. Anyway, it is obvious that the tails of VAF are very 
flat. Then the integration of VAF must produce a divergent result and T> becomes infinite, 
leading to the diffusion process anomalously enhanced. This result agrees with that obtained 
in the framework of the continuous-time random walk approach |19j (Levi walks) predicting 
the enhanced diffusion for power-law distributions of free paths; for (^) the motion becomes 
ballistic: T> diverges linearly with time. 
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Figure captions 

FIG. 1. The velocity variance calculated from Eq.@ (solid line) and resulted from both 
simulations: with the KP (dashed line) and using the deterministic random force (dots). 
The parameters: T = 1, m = 1 and e = 0.01. 

FIG. 2. The velocity autocorrelation function calculated using the KP with to = 1.5 
(dot-dashed line) and to — 3 (dashed line), normalized to unity at t — 0. The result of the 
simulation with the deterministic random force is marked by dots. The solid line shows the 
analytical result. The parameters are the same as in Fig.l. 
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